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ABSTRACT 



o 

^ ■ Context. New generations of instruments provide, or are about to provide, pan-chromatic images of debris discs and photometric measure- 
C/2 , ments, that requires new generations of models, to in particular account for their collisional activity. 

^ , ■ Aims. We present a new multi-annulus code for the study of collisionally evolving extended debris discs. We first wish to confirm and 
extend our early result obtained for a single-annulus system, namely that the size distribution in realistic debris discs always departs from the 
theoretical collisional "equilibrium" dN oc R^^ ^dR power law, especially in the crucial size range of observable particles (R ^ I cm), where it 
^ ' displays a characteristic wavy pattern. We also aim at studying how debris discs density distributions, scattered light luminosity profiles, and 
' SEDs are affected by the coupled effect of collisions and radial mixing due to radiation pressure affected small grains. 
Methods. The size distribution evolution is modeled over 10 orders of magnitude, going from yum-sized grains to 50km-sized bodies. The 
model takes into account the crucial influence of radiation pressure-affected small grains. We consider the collisional evolution of a fiducial, 
, idealized a=120 AU radius disc with an initial surface density oc a". Several key parameters are explored: surface density profile, system's 
\^ ' dynamical excitation, total dust mass, collision outcome prescriptions. 

Results. We show that the system's radial extension plays a crucial role and that the waviness of the size distribution is amplified by 
C inter-annuli interactions: in most regions the collisional and size evolution of the dust is imposed by small particles on eccentric or unbound 
orbits produced further inside the disc. Moreover, the spatial distribution of all grains S 1 cm significantly departs from the initial profile in 
^ E(a) oc a", while the bigger objects, containing most of the system's mass, still follow the initial distribution. This has consequences on the 
• ^ I scattered-light radial profiles which get significantly flatter, and we propose an empirical law to trace back the distribution of large unseen 
parent bodies from the observed profiles. We also show that the the waviness of the size distribution has a clear observable signature in the 
^ , far-infrared and at (sub-)millimeter wavelengths. This suggests a test of our collision model, that requires observations with future facilities 
. - - 1 such as Herschel, SOFIA, SCUBA-2 and ALMA. We finally provide empirical formulae for the collisional size distribution and collision 
timescale that can be used for future debris disc modeling. 
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1. Introduction 

Extrasolar discs around young stars have been imaged for more 
than two decades now. Observations and modeling have re- 
vealed the great diversity of these systems, in particular regard- 
ing the luminosity and density distribution of the dust com- 
ponent. Systems with the lowest dust to star l uminosity ratios 



have been commonly labeled debris discs (e.g. Lagrang e et al 



20001) . The archetypal member of this group is jSPictoris, 
which has been tho roughly observed and m odeled since the 
first observation by Smith & Terrile ('1984) (se e reviews by 



Vidal- Madj ar et al. , ,1994; .Kalas & Jewitt ,19951: lArtvmowicz 
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19971) . These systems are believed to represent a later stage 
of disc evolution, where most of the initial solid mass has al- 
ready been accumulated into planetary embryos or removed by 
collisional er osion and press u re forces (stellar ra diation/wind 



pressure, see GreavesI I2OO5I: iMeyer et al.1 bood for recent 
reviews on this subject). Simple order of magnitude esti- 
mates show that the dust in these system s cannot be primor - 
dial and has to be constantly replenished jArtvmowic 3 ll997h . 
Although cometary evaporation could also be a possibility 
(.Li & G reenberg 1998), the most likely dust produ ction mecha- 
nism i s collisional erosion of bi gger solid objects dArtvmowicz 
19971: iDominik & Decinll2003 i). This hypothesis is reinforced 



by the estimated ages of these syste ms, which are generally 
more than 10' yrs old ( Greavesll2005 ). For such ages, the stan- 
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dard planetary formation model ("e.g. lLissaueij|l993b predicts 
that most early stages of planetary formation, i.e. grain coag- 
ulation, planetesimal formation, runaway and/or oligarchic ac- 
cretion among these planetesimals, should already be over and 
that these systems should be made of large planetary embryos 
as well as smaller objects leftover from the formation process. 
The presence of big embryos should dynamically excite the 
system and lead to highly destruc tive mutual encounters be - 
tween the smaller leftover bodies dKenvon & Bromlevll2004l) . 
thus triggering a collisional cascade producing objects down to 
very small dust grains. 

The problems faced when modeling debris discs are numer- 
ous. One first difficulty is that all objects bigger than about 1 cm 
are completely undetectable by observations. Current observa- 
tions only probe the lower tail of a collisional cascade among 
objects invisible to us. The challenge is thus to reconstruct 
this hidden bigger object population from the observed dust 
component. But even for particles in the "observable" range, 
it is very difficult to get a coherent global picture. Each type 
of observations (visible, near-IR, far-IR, mm,etc...) is indeed 
predominantly sensitive to one particle size range and to one 
radial region of the disc. And even when a large set of such 
observational data at different wavelengths is available (includ- 
ing spatially resolved images, as for example for fi Pictoris), it 
does not allow to straightforwardly reconstruct the dust popu- 
lation. This "connecting the dots" procedure is always model 
dependent because it depends on many parameters, linked to 
the dust's composition, temperature, optical properties and size 
distribution, which can never be unambiguously constrained in 
a non-degenerated way (see for instan ce the thorough b est-fi t 



larger than 7?pr. The overabundance of Ri grains, in turn, in- 
duces a depletion of R2 objects with R2 slightly larger than 
R\, etc... This domino eff'ect propagates towards bigger sizes 
and leaves a characteristic wavy size distribution, with a pro- 



studies of Li & Greenberg ( 199% and lAugereau et al.l (120011) 
foryS Pictoris or lSu et al.l (I2005h for Vega). One challenge is in 
particular to get a coherent link between the mm-sized popu- 
lation, where most of the mass of the "dust" component is sup- 
posed to lie but for which spatial information is usually very 
poor, and the yum-sized grains, which should contain most of 
the optical depth and for which high-resolution observations 
are more and more frequently obtained. 

2. Previous works and paper overview 

The most basic way to perform these reconstructions of the un- 
seen big objects population or to derive coherent models of the 
dust population is to assum e that the classica l collisional equi- 
librium size distribution of iDohnanyi ( 19691) in dA^ oc R^^-^dR 
holds for all object sizes R. However, there are many reasons to 
believe that such an assumption is pro bably misleading. As it 
has been shown bv iThebault et al.l(l2003i, hereafter TAB03), the 
main problem arises from the smallest grains, whose behaviour 
is strongly affected by pressure forces imposed by the cen- 
tral star: radiation pressure in the case of lu minous stars, wind 
pressure for low-mass stars (e .g. AUMic, lAugereau & Beust 
2006HStrubbe & ChianjbOOeh . For stars of mass M, > 1 Mq, 
one major point is the presence of a minimum size cutoff 7?pr, 
all objects R < 7?pr being blown away by radiation pressure. 
Qualitatively, this depletion of R < /?pr grains leads to an over- 
density of slightly bigger grains Ri, because R < /?pr grains are 
depleted and can no longer efficiently destroy nor erode grains 



spect to the R power law (e.g. 


Cami 


po Bagatin et al. 1994; 


Thebault et al. 2003; Krivov et al. 


2006 


). These discrepancies 



with the dA^ oc R^^ ^dR distribution are reinforced by the fact 
that the smallest objects in the R > RpR range are put on very 
eccentric orbits by radiation pressure and have a dynamical be- 
haviour very different from that of the bigger no n radiation- 
pressure-affected bodies (see iThebault et alJ 2 



2003, for a thor- 



ough discussion on this topic). 

In TAB03 we quantitatively studied these complex effects 
for the specific case of the inner /3 Pictoris disc. For this pur- 
pose, a statistical numerical code was developed, which quan- 
titatively follows the size distribution evolution of a popu- 
lation of solid bodies, in a wide micron to kilometre size- 
range, taking into account the major effects induced by radi- 
ation pressure on the smallest grains (size cutoff, perturbed dy- 
namical behaviour,...). Our main result was to identify an im- 
portant departure from the dA^ oc R'^ ^dR law, especially in 
the 1 fim to 1 cm range. The main limitation of this code is 
that it considers a single, isolated annulus. It can thus only 
be used to study a limited region at one given distance from 
the star (^ 5 AU in the case considered) but not the system as 
a wh ole. A multi-an nulus approach is needed to achieve this 
goal. iKenyon & Bro mley (2002, 2004) have developed such 
statistical multi-annulus codes, which have been applied to var- 
ious contexts. These codes are in some respect more sophisti- 
cated than the one used in TAB03, in particular because they 
follow the dynamical evolution of the system (which is fixed 
in TAB03). Nevertheless, the price to pay for following the dy- 
namics is that the modeling of the small grain population is 
very simplified, with all bodies below a size R ^ Im following 
an imposed R 'ldR size distribution, thus implicitly overlooking 
the aforementioned consequences of the s pecific behav i our of 
the smallest dust particles. More recently. iKrivov et al.l (l2006h 
used a different approach based on the kinetic method of sta- 
tistical physics. This model is able to follow the evolution of 
both physical size and spatial distribution (ID) of a collision- 
ally evolving idealized debris disc, from planetesimals down 
to yum-sized grains. This model has also the added advantage 
of taking into account a large range a unbound particles below 
the blow-out limit. This innovative approach gave promising 
results for the specific case of the Vega system. However, the 
modeling of collisional outcomes is, as acknowledged by the 
authors themselves, very simplified, with for instance all cra- 
tering impacts being neglected. 

In this paper we present a newly developed multi-annulus 
version of our code, aimed at studying the collisional evolu- 
tion of spatially extended systems. Intra and inter-annuli in- 
teractions, due to the radial excursions of radiation-pressure 
affected small grains, are considered. In addition to this new 
global scheme, a new and improved modeling of collision out- 
comes is presented, with particular attention paid to the crucial 
cratering regime (Section |3] and Appendix). In order to clearly 
identify and study the complex mechanisms at play, we con- 
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sider in the present study the case of a fiducial ideahzed de- 
bris disc of 120AU radial extension, and explore surface den- 
sity distributions in oc a" around the reference MMSN 
a = -1.5 case, where a is the distance to the star The evo- 
lution of the system's size distribution, and its significant de- 
parture from the standard Dohnanyi steady-state power law, 
is followed until t = lO^yrs and is presented in section |4] 
The role of several key free parameters, such as the system's 
dynamical state, stellar mass and grain physical composition 
are explored in section |5] The evolution of the system's spa- 
tial distribution, optical depth and the correspondence between 
observed dust and unseen bigger parent bodies is addressed in 
section |6] In section [7] we investigate the impact these results 
have on important observables, in particular the scattered light 
and thermal emission luminosity profiles as well as the SEDs. 
In section[8]we discuss the robustness of our results and derive 
empirical laws for the size distribution and collisional lifetimes 
which might be extrapolated to any kind of extended collision- 
ally evolving debris disc. Conclusions and perspectives are pre- 
sented in section|9] More specific studies of specific debris disc 
systems will be the purpose of a forthcoming paper 

3. Numerical model 

3.1. Structure 

Our code adopts the classical particle-in-a-box statistical ap- 
proach to follow the collisional evolution of a population of 
solid bodies of sizes comprised between /?min < ^pr, where 
RpR is the radiation pressure blow-out size, and /?n,ax - 10- 
50 km. The system is made of concentric annuli of width 
Afl,fl and centered at distances «,„ from the star Within each 
annulus ia, bodies are distributed into n size bins, each bin cor- 
responding to bodies of equal size Ri. The evolution of the size 
distribution with time is given by the estimated collision rates 
and outcomes between all coUisionally interacting (ia, i) bins. 
For small particles produced in an annulus ia and placed on 
high-eccentricity or unbound orbits by radiation pressure, col- 
lisions with bodies located within all annuli crossed by their or- 
bits are taken into account. A detailed presentation of the model 
is given in AppendixlAl 

One key parameter for our model, and any similar study 
for that matter, is the prescription for the collision outcomes. 
We adopt the classical approach where the outcome of an im- 
pact between a target of size R, and a projectile of size Rp de- 
pends on the ratio between the center of mass specific kinetic 
energy of the colliding bodies E^oi and the the so-called critical 
specific shattering energy Q*, which depends on the objects' 
sizes and composition. Depending on the respective values of 
these two parameters, impacts result in catastrophic fragmenta- 
tion, cratering or accretion. The collision-outcome prescription 
has been updated with respect to the one in TAB03, in partic- 
ular for what concerns the cratering regime. The new model 
now also accounts for differential chemical composition within 
the system, the main parameter being here the radial distance 
from the star Ojce below which water ice sublimates. The com- 
plete collision outcome procedure is described in more details 
in AppendixlB] As explained in this Appendix, we consider a 



"nominal" case for the fragmentation and cratering prescrip- 
tions and with ai^s = 20 AU, but other cases are explored (see 
section|531i. 

The price to pay for following the size distribution over 
more than 10 orders of magnitude in size is that we can- 
not accurately follow the dynamical evolution of the sys- 
tem, whose dynamical characteristics have to be fixed as in- 
puts. In this case, all CPU-time consuming calculations of 
mutual impacting velocities and collision physical outcomes 
are performed onc e at the beginning of the run (e.g. TAB03, 
iKrivov et al. 20061) . We shall therefore implicitly assume that 
the disc has reached a quasi-steady dynamical state, which 
holds for timescales longer than the ones considered in the sim- 
ulations. We consider identical average values for particle ec- 
centricities and inclinations for all size bins, with the exception 
of bins corresponding to particles affected by radiation pres- 
sure for which specific orbital characteristics are numerically 
derived (see Appendix). 

For a more detailed description of our code, see the 
AppendixieslAl and |B] 

3.2. Initial conditions 

As mentioned in previous sections, we consider here a fiducial 
idealized debris disc, for which the initial spatial distribution 
follows the classi cal Minimum M ass Solar Nebulae (MMSN) 
profile derived by Hayashij (Il98lh . where the surface number 
density is such that £(«) oc where a is the distance from 
the star We consider all concentric annulus disc, that extends 
from flmin = 10 AU to a,nax = 120AU, a typical range for the 
radial extension of dusty debris discs. 

The initial conditions are chosen in accordance with the 
current understanding of debris discs, i.e. systems in which 
the bulk of planetesimal accretion process is already over and 
large planetary embryos are present. These large objects should 
dynamically excite the system, and average eccentricities and 
inclinations in the disc ma y reach values of t he order of 0.1 
for Lunar-sized embryos (lArtymowic j 1 1 997h . We thus take 
(e) = 0.1 = 2{i} (with (/) in radians) as our nominal dynam- 
ical conditions and explore diff'erent orbital values in separate 
runs. We follow the collisional evolution of all objects in the 
lRmm,Rm!ix] range, where R^i„ < /?pr and R^,^^ ^ 50km. We 
take as a reference value Rpr = 5 fim, which corresponds to the 
value for a compact grain around a ySPictoris -like star (A5V), 
but other possible 7?pr values for earlier and later type stars are 
also explored (section lS^ . The planetary embryos themselves 
are left out of our study since they are too isolated to contribute 
to the continuous collisional cascade, and can only affect the 
dust production rate through sudden is olated events (for the de - 
tailed study of such violent events, see lGrigorieva et all 2007 ). 
We shall assume that the initial size distribution at t = Oyr 
follows the classical dA^ oc R^^ ^dR power law from R^in to 
^max and we follow subsequent departures from this "equilib- 
rium" distribution as time goes by. Having fixed the initial 



' However, the initial size-distribution is not a crucial parameter, 
since test runs have shown that the system always settles to the same 
steady state regardless of the initial dN/dR prescription. 
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Table 1. Nominal case setup. The fields marked by a y are 
explored as free parameters in the simulations. See text for de- 
tails. 





Radial extension 


10 < a < 120 AU 




Number of annuli x radial width 


11 X lOAU 


V 


Initial surface density profile 


E(a) oc a"'-^ 


V 


Total "dust" mass (R < 1 cm) 


0.1 Me 




Size range modelled 


3 fj^m< R < SOkm 




Number of size bins 


103 




Initial size distribution 


dA' oc R-^ ^dR 


V 


Sublimation distance (water ice) 


flsiib = 20 AU 


V 


Dynamical excitation 


<e) =0.1 =2(0 


V 


Stellar type 


A5V 


V 


Blow out size 


RpR = 5 fim 


V 


Collision outcome prescription 


(see Appendix int 



size-distribution, the initial disc mass is a free parameter which 
is explored in separate runs. This disc mass is parameterized 
by Mdust, the total amount of "dust", i.e. grains smaller than 
^ 1 cm, in the system. The parameter Mjust has been chosen as 
a reference because it is usually the most reliable constraint on 
the disc mass which can be derived from observations, since 
larger objects are observationally undetectable. Most of this 
dust mass is believed to be contained in the bigger millimetre- 
sized grains detected at sub-millimetre to millimetre wave- 
lengths. Such millimetre wavelength surveys have shown that 
for debris discs around young main sequence stars, Mjust is 
typicall y comprised be tween 0.001 and a few 0.1 M© (e.g. re- 



view by lGreavesll2005L and references therein). Accordingly, 
we shall consider two limiting cases: a low mass disc with 
^du.st = O.OOlMe, and a high mass disc with Mdust = O.IM© 
(in both cases, the initial distribution of bigger objects is ob- 
tained by extrapolating a dA^ oc R^^ ^dR size distribution up to 
^max)- Particles within the size bins are assumed to be compact 
silicates in the regions closer to the star than the subimation 
limit Osub and compact ices beyond Osub, with Usub - 20AU in 
the nominal case (see section B.l of the Appendix). 

For each run, we let the system evolve for 10^ yrs. Of course 
debris discs can have ages exceeding by far this value (as for 
instance Vega or e Eridani), and longer timescales should in 
principle be considered here. We should however restrict our- 
selves to 10^ yrs because in most of the cases the system reaches 
a steady-state much earlier than this (typically after ~ lO^years 
for our nominal case). The only exception to this is the "low- 
mass" case with Mdust - 0.00 IM©, for which the steady state 
is not reached, at ffinai = lO^yrs, in the outer regions of the 
systems. For this specific case, we let the system evolve until 
lO'^yrs. 

All initial parameters for the nominal high mass case are 
summarized in Table[T] 

4. Results for the nominal case 

4.1. High-mass disc fMdust = 0.1 M^j 

Figure [T] shows the evolution of the size distribution for four 
annuli at different distances from the star. In the innermost an- 



nulus (Fig.[T^), a weak wavy pattern develops, starting with the 
depletion of R < /?pr grains and propagating upward. Once 
the pattern has fully developed, subsequent evolution consists 
in a progressive total mass loss while the global size distribu- 
tion profile is conserved. This wavy pattern is however much 
less pronounced in this innermost annulus than in TAB03. The 
first reason for this is that TAB03 considered a region further 
inside the disc, at 5 AU, while the present inner annulus starts 
at a = 10 AU and extends up to 20 AU. Impact velocities, and 
their destructive efficiency, are thus significantly lower here. 
The second reason is due to our revised collision-outcome pre- 
scription, in particular for crateiing events, which in TAB03 
had a dominant role in shaping the size distribution in the 
R $ I cm domain (see Table 4 of this paper). With the more 
realistic cratering prescription taken here, excavated masses 
are significantly smaller in the small grains domain than in 
TAB03 (see Sec. IB. 3b . hence the shallower patterns in the size 
distribution. The knee in the distribution around 0.1-1 km is a 



well known feature (e.g Campo Bagatin et al.lll994l) due to the 
switch from the strength dominated regime, where bodies re- 
sistance weakly decreases with increasing size, to the gravity 
dominated regime, where bodies resistance to impacts rapidly 
increases with increasing size. It can be easily checked that the 
location of the knee at R ^ 0.1km coiTesponds to the least 
impact-resistant bodies (see Equ. lB.lb . Furthermore, for large 
objects, reaccumulation of fragments after an impact also be- 
gins to play a major role. 

In the more distant annuli, on the contrary, very pronounced 
wavy patterns are observed in the size-distribution (Figs.[TJ),c 
and d). The most striking features are the overdensity of R ^ 
1 .5 7?pR bodies, and above all, the strong depletion of bodies in 
the submillimetre range (^ 10-50 7?pr). This result might ap- 
pear counter-intuitive since one would expect these features to 
be even less pronounced than in the innermost annulus because 
of the longer dynamical timescales and lower impact veloci- 
ties in the outer regions. The main cause for these sharp fea- 
tures are in fact small high-yS grains originating from other an- 
nuli further inside the disc (where /3 classically designates the 
radiation pressure to gravitation forces ratio). This is clearly 
illustrated in Fig. [2] which compares, in the middle 50-60 AU 
annulus, the final size distribution (solid line) with the size 
distribution obtained when only considering locally produced 
particles (dashed line). In the R ^ 30/im range, foreign-born 
grains make up up to 90% of the local population, thus result- 
ing in a factor ^ 10 increase of the number density. But the 
effect of these additional inner-disc-produced grains on the 
system's evolution exceeds by far that simply due to a num- 
ber density increase of an order of magnitude. Indeed, as these 
grains have had more time to reach high radial velocities than 
the locally produced grains of the same size, they will impact 
objects in the annulus at much higher relative velocities. As an 
example, for a target on a circular orbit at 50 AU, an impact 
by a locally produced small grain with /5 - 0.45, will occur at 
Av ^ 1 km.s"', whereas an impact hy a /3 - 0.45 grain pro- 
duced at 10 AU will occur at Av ^ 5 km.s This will result in 
much more destructive collisions. It is this higher destructive 
efficiency which is responsible for the deep depletion of objects 
up to - IQQRpR - 0.5 mm. Another important result is that a 
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jg30 I , 1 , ^ , 1 , 1 , ^ , 




size (cm) size (cm) 



Fig. 1. High-mass case (Mdust - O.IM^). Time-evolution of the size distribution at four different locations in the disc. Note that 
the y-axis displays the mass contained in one size bin, which is a correct way of displaying the mass distribution since all size 
bins are equally spaced in a logarithmic scale. 



large fraction of the sub-mm grain depletion is due to cratering 
impacts, as appears clearly from the test run with no-cratering 
shown in Fig. |2] (dotted line). Indeed, small R < 30fim grains 
cannot directly break-up objects bigger than ^0.1 mm, even 
for their increased impact velocities, while they can efficiently 
erode by cratering bodies up to almost ^ 1 cm. 

4.2. Low-mass disc (M^ust = 0.001 M®J 

The evolution of the size distribution for the low-mass case is 
displayed in Figure |3] The main difference with the Mdust = 
0.1 Me case is that the global evolution of the system is much 
slower The slowing down is logically of the order of the discs 
mass ratio (i.e. a factor of about 100). In the a ^ 70 AU region, 
after 10^ years, the system has reached a quasi steady-state rel- 
atively similar to the high mass case, with an overdensity of 
^ 1 .5 RpR grains, followed by a depletion of sub-mm grains of 
approximately one order of magnitude compared to the initial 
size distribution. 



In the outermost regions, however, we observe a much 
deeper depletion of sub-mm grains than in the high-mass case. 
This is because the collisional-equilibrium, contrary to the in- 
ner disc regions, has not been reached at f = 10^ years: the 
erosion of sub-mm grains, by high-/? particles coming from 
the inner regions, has already reached full efficiency (after less 
than lO^yrs), while the production of new sub-mm grains by 
erosion of larger objects occurs on much longer time-scales, 
exceeding lO^yrs in the outer regions (see more detailed dis- 
cussion in the next section). This is clearly illustrated in Fig.[3}l, 
which shows that, in the outermost annulus, the population of 
> 1 cm objects remains largely unaffected by collisional pro- 
cesses after 10^ years. In order to get a better idea (and despite 
of the huge CPU-time cost), we decided to let this low-mass 
disc collisionally evolve for another f - 10** years. As can be 
seen in FiglJl, at this later time the quasi-steady state is al- 
most reached in the outermost annulus, but the second "knee" 
in the size distribution atR ^ 0.1km is still not visible yet. The 
full steady-state is here probably reached on timescales of the 
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'locally produced grains only 



0.0001 0.0010 0.0100 0.1000 1.0000 10.0000 100.0000 

size (cm) 

Fig. 2. Size distribution for the 50-60AU annulus, at f = 
10^ yrs, for the nominal case (solid line), when only taking into 
account fragmenting impacts, i.e. no cratering (dotted line), and 
when only taking into account the locally produced grains, i.e. 
no impact with grains coming from inner annuli (dashed Une) 



order of ~ 1 Gyr, which are presently out of the reach of our 
numerical code. 



4.3. Collisional lifetimes 

We define the collisional lifetime of a particle as the average 
time fcoii it takes for the object to lose 100% of its mass by col- 
lisional processes. Let us point out that the collisional mass loss 
has two origins: (i) catastrophic fragmentation, for which a par- 
ticle loses by definition 100% of its mass at each fragmenting 
encounter, and (ii) cratering, for which the particle is progres- 
sively eroded after each impact (excavated mass M„a given in 
App. lB.3l l. The left and right panels of Fig. |4] display the values 
of fcoii(a,/?), after 10^ years, at different radial locations in the 
disc, for both the nominal high-mass and the low-mass cases. 
Note that for particles placed on high-eccentricity orbits by ra- 
diation pressure, fcoii(fl,/?) is the collisional lifetime of a parti- 
cle initially produced at distance a, when taking into account 
all the collisions this particle will suff'er in the different annuli 
it will cross on its eccentric orbit. 
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Fig. 4. Collisional lifetimes of the particles as a function of their size, at f = 10^ years, and at four different locations in the disc: 
a = 15 AU (solid line), a -25 AU (dotted line), a = 55 AU (dashed line), a = 115 AU (dot-dashed line). The four horizontal lines 
are, for each distance, the collision timescales deduced from the simplified T^^^^ia) - (tQ)"' formula. The a values refer to the 
center of the annulus where the particle has been produced. Left panel: high-mass case (Mdust = 0.1 Me), right panel: low-mass 
case (Mdust = 0.001 Me) 



We discuss first the case of the high-mass disc (Fig.lH left 
panel). For large objects, we obtain the predictable result that 
collision lifetimes increase with increasing distances from the 
star. This results from three concurring factors: particles num- 
ber densities decrease with a, dynamical timescales get longer, 
and impact velocities lower, leading to less eroding impacts. 
For objects in the dust size range, however, the situation is 
much more complex, mainly because of the major influence 
the radial movements of small high-/? grains have on the col- 
lisional evolution. We find that the most short-lived particles 
are the ones with R ^ 100//m, which logically corresponds to 
the most depleted population in the system (Fig.[TJ. For grains 
with sizes R < 100/im, fcoii very rapidly increases with de- 
creasing sizes. The explanation for this trend is twofold: first, 
destructive impactors in the R < Rp^ range are strongly de- 
pleted, and second, many of these high-/? grains spend a large 
fraction of their eccentric orbits in the empty region of the disc 
beyond 120 AU. This global trend of the fcoU dependence with 
size could to some extent be compared to the one obtained by 
[Stmbbe & Chiang (2006) for the AU Mic system (where stellar 
wind from the central M-type star could pl ay the same role as 
radiation pressure around A-type stars, see lAugereau & Beusli 
[2OO6). These authors also found dfcou/d/? > for large grains 
and a very sharp dfcou/d^ < gradient for smaller grains (see 
Fig.l of this paper). However, these similarities are only qual- 
itative (with major quantitative diff'erences regarding the turn- 
off size or the slope of the dfcnii /d R laws) and should in any c ase 
be taken with great care, since the lStrubbe & Chiang! ( 200 6^ es- 
timates were obtained for a radially narrow system and with a 
simplified analytical law for the collision rates and outcomes. 



is here no flux of destructive small (high-yS) impactors com- 
ing from further inside the system, contrary to the other an- 
nuli. Note that the only objects having fcoU > ffinai = lO^yrs 
are the largest planetesimals, of size R > 0.2 km in the outer- 
most regions, and R > 3 km in the inner annuli. This means 
that, as a first approximation, all sub-kilometre sized objects 
are collisionally evolved, i.e., no object other than the largest 
kilometre-sized bodies are primordial. 

All these global trends are also valid for the low-mass case 
(Fig.m right panel). However, the fraction of primordial ob- 
jects is much higher than in the high-mass run. In the outer 
annulus, for example, no object bigger than ^ 1 cm has been 
collisionally processed in 10^yrs|3. This is in agreement with 
what was pointed-out in the previous section, namely that the 
collisional cascade did not fully develop in the outermost re- 
gions after lOMyr of evolution. For the dust-size range, how- 
ever, the result that all particles are collisionally processed dur- 
ing the system's lifetime still holds. This is why the shape of the 
size distribution is relatively similar in the high and low-mass 
runs. 

As an interesting comparison, we also plotted on the graphs 
the collisional timescale obtained by the formula f*'^jj - (rfl)"', 
where r is the geometrical vertical optical depth and Q the an- 
gular velocity. This simplified relation is indeed often used in 
the literature as giving an approximate estimate of the colli- 
sion lifetimes of the smallest grains. As can be clearly seen, it 
proves to be a very poor match to our numerically derived life- 
times. Diff'erences can reach up to 2 orders of magnitudes in 
the crucial R < 0.1 mm range. 



The relative lifetimes between different regions of the disc 
follow the logical trend in dfcou/da > 0, with the important 
exception of the innermost annulus, for which collisional life- 
times of small grains are relatively high, simply because there 



- For sake of comparison with the high-mass case, we do here con- 
sider the same fy,„„; = 10'' years value for the low-mass run, instead of 
the additional lO^years explored in Fig.3 
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Fig. 5. Impact of the disc's dynamical excitation on the size 
distribution after 10^ years of evolution. The figure shows the 
size distribution of the whole system (i.e. all 1 1 annuli) assum- 
ing an initial dust disc mass of Mdust = 0.1 M® (high-mass 
case). 



5. Parameter dependence exploration 

5.1. Dynamical excitation: (e) 

The exact orbital distribution of particles in debris discs is in 
general very poorly constrained. The only observational con- 
straint comes from measuring the disc's vertical thickness and 
deriving estimates of orbital inclinations, but such constraints 
are scarce. Edge-on discs represent the most favorable cases 
since H/a, where H denotes the vertical scale height, can be 
directly measured. Five out of a dozen of spatially resolved 
discs have this particular orientation: P Pictoris , AU Mic, 



HP 3 2297 (ISchneideretal.1 l2005i) HD 139664 (,Kalas et al 



20061) . HD 15115 dKalas et al.ll2007h . For th e two most s tudied 



discs, only partial information is available. iKrist et al.l (I2005h 
find H/a ^ 0.04 in the case of the AUMic disc, with ratios as 
small as 0.02 close to the position of maximum surface density. 
The /3 Pictoris disc appears geometrically t hicker with H/a ra- 
tios as large as - 0.1 (.Golimowski et al ]|2006) . However, these 
meas urements include the so-c alled disc warp which, accord- 
ing to lGolimowski et al. I (l2006h . might be due to a blend of two 



separate, intrinsically thinner disc components inclined with re- 
spect to each other by a few degrees. The Pictoris disc might 
then in fact be less vertically extended than it appears to be. 
The modeling and inversion of scattered light brightness pro- 
files of inclined, ring-shaped discs do not provide much more 
constraints. The HR4796 and HD 181327 rings for example, 
might have H/a ratios as large as about 0.1 at the positions of 
maximum sur face density, but the actual ratios could be two 
times smaller ( Augereau et alj[l999l : Schneider et al llaooe*). As 
pointed out in section [331 other estimates of the disc's verti- 
cal thickness come from general theoretical arguments. Debris 
discs are indeed believed to correspond to the late stages of 
planetary formation where Lunar-to-Mars sized embryos dy- 
namically excite the system. However, this argument can only 



lead to rough order of magnitude estimates of the dust's orbital 
elements. It is thus important to explore different possible val- 
ues of (e) and (/). Due to the CPU-time consuming aspect of 
the simulations, we chose to restrict ourselves to the high-mass 
system and perform one additional "dynamically colder" case 
with (e) - 0.03 = 2(/), one "very cold" system with (e) - 0.01 
and one dynamically "hotter" case with (e) = 0.2. A compari- 
son between these three cases and the nominal (e) = 0.1 case 
is displayed in Fig.|5] For sake of clarity, we consider here the 
whole system, summing up the contributions of all radial an- 
nuli. 

Contrary to what could be intuitively expected, the deple- 
tion of objects smaller than 1 mm is more pronounced in the dy- 
namically cold case (e) = 0.03 (Fig.|5]l. There are two concur- 
ring explanations for this apparent paradox. On the one hand, 
the rate at which sub-millimetre grains are eroded only weakly 
depends on the system's dynamical excitation. Indeed, the ve- 
locity at which these grains are impacted by smaller micron- 
sized particles is mainly imposed by the strong radiation force 
acting on the latter and only weakly depends on the eccentric- 
ity of their parent bodies' orbits. On the other hand, the rate at 
which big grains are produced, by impacts between larger ob- 
jects, strongly depends on the system's dynamical excitation, 
since these larger objects' orbits, and thus their impact veloc- 
ities, are insensitive to radiation pressure effects. As a conse- 
quence, the balance between production and erosion of sub- 
mm grains is more negative for low (e) values of the parent 
bodies orbits, hence the more pronounced depletion. For the 
"very" cold case, this effect is even more pronounced, and one 
can witness a global general depletion of all dust size grains, 
while objects in the > 1cm range are mostly unaffected by any 
collisional evolution. 

For the high-excitation case, the depletion of sub-mm 
grains is almost identical to the nominal case, which is here 
again a direct consequence of the fact that the dynamics of the 
very small grains is controlled by the radiation pressure force. 
These results are clearly illustrated in Fig.|6] showing the colli- 
sional lifetimes in both the dynamically "hot" and "cold" cases. 
While fcoii is roughtly inversely proportional to (e) for large 
(^ 1cm) particles, the collisional lifetimes of small grains only 
weakly vary with the average dynamical excitation in the disc. 

5.2. Mass of tine star and Rpr value 

The nominal case considered in our simulations is that of a /3- 
Pictoris like star of mass M* = 1.7 M© and a corresponding 
radiation pressure cut-off size RpR = 5 fj.m. We explore here 
the M, and /?pr parameters by considering, in addition to the 
nominal case, two limiting cases: one G-star of mass 1.1 Mq 
with /?pR = 1 jim, i.e., the lowest star mass for which compact 
silicate grains can reach the /3 = 0.5 limit, and one Vega-like 
AOV star of mass 2.5 M© and Rpk - lOAtm. All /^pr values 
have been derived using the lGrigorieva et aL (2007) algorithm. 

As appears clearly on Fig.|7] the size distributions for all 
three systems in the "dust" grains size range (R < I cm) are 
relatively similar. The profiles are shifted in size with respect to 
each other, reflecting the difference in i?pR values. Interestingly, 
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Fig. 7. Impact of the radiation pressure cut-off size value 7?pr 
on the size distribution of the whole system (i.e. all 1 1 annuli) 
at f = 10^ years for the high-mass case (Mjust = 0.1 Me). 
Solid line: nomical case, M, = 1.7 Mq and 7?pr = 5 yum, dotted 
line:M^ = 1.1 and 7?pr = 1 jim, dashed line: M* = 2.5 Mq 
and PR = lOyum 



the location of the overdensity of smallest grains is always 
given by the relation R ^ 1.5 RpR, while the most pronounced 
depletion is always obtained for R ^ 10-50 RpR. However, the 
amplitude of this depletion increases with increasing RpR val- 
ues (i.e. star masses). This is because, in the strength regime, 
smaller grains are more resistant to impacts than bigger ones, 
which implies that an impact between, say, aR - 1 .5 RpR dust 
grain and aR - 100 RpR object is more erosive for larger values 
of 7?pR. Furthermore, for a more massive star, impact velocities 
are higher (for the same orbital parameters), which also leads 
to more destructive collisions. 
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Fig. 8. Size distribution, at f = 10^ years, for the whole system 
(all 1 1 annuli), for 4 different initial surface density distribu- 
tions 2(a) oc a", keeping the system's total mass constant. 



5.3. /n/f/a/ density profile 

We have considered as a standard case a system following a 
standard MMSN spatial distribution in 'E{a) oc o^'^. However, 
in order to check the robustness of our results, other indexes 
for the E(fl) oc a" dependence have been explored. Fig.[8]shows 
that the global size distributions within the system only weakly 
depends on the initial 2(fl) power law. The only noticeable 
trend is a slight damping of the wavy distribution in the R 1cm 
range for flatter S(a) profiles. This result is logical since we 
have seen in section 14.11 that, in a given region of the disc, 
the evolution of the sub-mm grains is mainly imposed by the 
flux of high-yS particles coming at high radial velocities from 
the inner regions. The influence of these inner-disc born grains 
should logically diminish for less steep l,{a) profiles, for which 
their relative abundance compared to the local population is 
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smaller. However, these differences between the different 2(a) 
cases remain limited in amplitude and all size distributions re- 
main very close to the result of the nominal case. 

5.4. Collision outcome prescription 

As discussed at length in the AppendixiBl the collision out- 
come prescription is a poorly constrained parameter, first be- 
cause of uncertainties regarding the chemical composition and 
structure of the grains and planetesimals in debris discs, and 
second because of significant differences between the predic- 
tions of all existing models. Our nominal case assumes a sub- 
limatio n distance for ices cisub = 20 AU, the Benz & Asphaug, 
(Il999h prescription for the critic al specific energy Q* for sili- 
cates and 2*,^ = 2sii/5, and the iKoschny & Griinl (1200 ih for- 
mula for crater-excavated masses M^ra for ices and silicates 
(see AppendixlBli. In order to explore how our results depend 
on the collision prescription, we have performed the two fol- 
lowing additional runs: 



- one "h ard " material run , assuming the iBenz & Asphaug 
(Il999h and lKoschnv & G riin (2001J) prescriptions for com- 
pact silicates hold for the entire disc (i.e. numerically set- 
ting flsub = AU) 

- One "weak " material run, wher e we assume the Q* pre- 



scription of iKrivov et al.l (l2006h . and a value of five 



times higher than in the nominal case. 

The results are displayed in Fig.|9] 

As could be logically expected, the wave-like structure is 
much less pronounced for the "hard" material run. As a mat- 
ter of fact, only the first wavy feature, affecting the smallest 
grains, is clearly visible, and its amplitude is damped by a 
factor ^ 3 compared to the nominal case. Moreover, the size 
for which the strongest depletion is reached is shifted from 
R - 20RpR = 1 00 yum to 7? ^ 4Rpr = 20yum. For the "weak" 
material run, the exact opposite is observed: pronounced wavy- 
features propagate up to the largest sizes, and the amplitude 
of the depletion of sub-mm grains is significantly increased 
and reaches almost two orders of magnitude. Contrary to the 
hard-material run, the depletion is now shifted towards bigger 
grains as compared to the nomin al run. The w e ak-m aterial run 
partially resembles the results of Krivov et al. ( 20061) . which is 
logical considering that we took identical Q* values, but differ- 
ences are observed, which can probably be attributed to the fact 
that cratering impacts are here taken into account. 

A comparison between Fig.|9] and all other parameter ex- 
ploration runs of Figs.|5] to [8] clearly shows that the collision- 
outcome prescription is the most crucial parameter the final 
size-distribution depends on. Unfortunately, this parameter is 
probably the most poorly constrained in the present problem. 
As described at length in the Appendix, particular attention 
has been paid here to this crucial issue. We have tried to im- 
prove on most previous studies (including TAB03) and con- 
sider an upgraded model incorporating the most relevant avail- 
able data for the Q* as well as fragmentation and, more specif- 
ically, cratering prescriptions. Nevertheless, large uncertainties 
remain. Firstly, important grain properties, which are crucial 
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Fig. 9. Impact of the outcome collision prescription on the size 
distribution at f = 10^ years and for the whole system (i.e. all 1 1 
annuli). Three different collision-outcomes prescriptions have 
been assumed: nominal case (solid line), "hard material" case 
(dashed fine) and "weak material" case (dotted line). See text 
for details. 

for understanding their response to impacts (ice fraction, poros- 
ity, etc.), remain poorly constrained for most debris discs. 
Secondly, even if all grain characteristics were fully known, 
it remains to see to which extent collision outcome energy- 
scaling models (even the more advanced version considered 
here), mostly obtained by experiments on cm-to-decimetre 
sized targets, might apply over such a wide size range, es- 
pecially for very small micron-sized grains. There is to our 
knowledge no fully reliable data on what the outcome of a colli- 
sion between, say, a 5fim grain and a 0. 1mm target at 500m. s"' 
"really" is. Basically, it all comes down to how soft or hard 
(with respect to a collisional event) particles in the < 1cm 
range are, and how these characteristics might vary with size. 
In this respect, we believe our nominal case collision prescrip- 
tion to be the most reliable one given the (still limited) cur- 
rent knowledge on this complex problem. Nevertheless, sig- 
nificantly different collisional behaviours cannot be ruled out. 
Fig. Improbably gives a good idea of realistic boundaries for the 
limiting "hardest" and "weakest" material cases, showing that 
the waviness of the size distribution decreases with increasing 
collisional resistance of the objects. 

6. Spatial distribution and dust to planetesimals 
mass ratios 

6.1. Radial distribution 

For sake of clarity, we consider here only the nominal high- 
mass run. FiglTO]clearly shows that the spatial distribution sig- 
nificantly departs from the MMSN profile for all objects in the 
"dust" size range (< 1cm). As could be logically expected, 
the strongest departure from the initial MMSN profile is ob- 
tained for grains in the sub-mm size range. For this population, 
the sharpest feature is a density drop in the regions just out- 
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Fig. 10. High-mass disc. Radial distribution, at"r=^ LOjyea 
of the mass surface density for different object sizes. For each 
size range, all surface densities are renormalized to the surface 
density in the first annulus. The dashed line represents the the- 
oretical distribution should a MMSN power law in a hold 
starting at the innermost annulus. 
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Fig. 11. High-mass disc. Radial distribution, at f = 10^ years, 
of the geometrical vertical optical depth for different size 
ranges, parametrized by their p parameter. 



side the first annulus. This drop is easily understandable and is 
due to the inter-annuli interactions already described in 3.1.1: 
in the innermost annulus, only locally produced small grains 
can erode sub-mm particles, but such locally produced small 
grains, blown out by radiation pressure on unbound or very el- 
liptical orbits, have not the time to be accelerated to high ve- 
locities, which limits their destructive or erosive power In all 
other annuli, on the contrary, small grains coming from the in- 
ner regions impact local bigger grains at very high velocities 
and are able to deplete them more significantly. 

For small grains in the < 50;um range, the radial distri- 
bution is very flat, even flatter than the one which should be 
expected in a steady flow of outgoing unbound particles, where 
simp le mass conser vation considerations lead to Z(fl) oc a"' 
(e.g. lSu et al.ll2005h . This profile cannot be explained by sim- 
ple blow out of unbound particles since most of the grains in 
the < 50 im\ range are on bound orbits {Ry-r - 5 fiva for our 
nominal case). On the other hand, the mass surface density dis- 
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Fig. 12. High-mass disc. Evolution of the vertical optical depth 
profile with time. 
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tribution of the total system (all particle sizes) is still relatively 
close to a classical MMSN profile in a"'^ (solid black line in 
Fig. [Toll. This is not an unexpected result, since the bulk of the 
disc's mass is still contained in the biggest, kilometre-sized 
particles, which are only marginally affected by specific col- 
lisional behaviour of the smallest grains. Therefore, there ex- 
ists a major discrepancy between the spatial distribution of the 
largest undetectable objects and that of the grains in the dust- 
size range, i.e. those accessible to observations. 

Another interesting result concerns the geometrical vertical 
optical depth T{a). Fig.[TT]shows the respective weight of differ- 
ent grain populations. We see that, except for the innermost re- 
gions, T(a) is completely dominated by grains from a very nar- 
row size range of a-meteoroids just above the blow-out limit 
y6 = 0.5 o 7? = /?pR. Of course, even with a standard power 
law distribution in dA^ oc R^^ ^dR, the optical depth should be 
dominated by small objects, since 

However, this tendency is much more pronounced here. As a 
matter of fact, when averaged over the whole system, it can be 
shown that 50% of the total optical depth is due to bodies in 
the 0.3 < j8 < 0.5 o R^r < R < 1.6Rpr range. For a Dohnanyi 
profile, the size range containing 50% of the total optical depth 
is much broader: RpR < R < 4/?pr. The temporal evolution 
of the T(fl) profile is also of interest. As Fig. [12] clearly shows, 
it rapidly settles (in a few lO^yrs) to a relatively "flat" radial 
profile, much flatter than the initial a one. This flattening 
is due to several mutually connected factors. The main one is 
due to what has been previously outlined, namely that the opti- 
cal depth is dominated by grains from a narrow size range just 
above Rp^. These very small grains are very quickly placed on 
very eccentric orbits, and will thus spend most of their orbits 
outside their annulus of production. As a consequence, small 
high-/? grains will naturally tend to be depleted in the inner re- 
gions and pile-up in the outer ones. In addition to this, the col- 
lisional erosion of bigger dust grains in the ~ 0.05 mm to 1 mm 
range, which make up most of the mass "reservoir" from which 
smaller high-/? grains are collisionnaly produced, is faster in the 
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Table 2. Relative mass fraction M^,„ contained in the smallest 
(R < 20 fim) grains, M„„„ in the biggest dust particles (0. 1 mm< 
R < 1 cm), and Mhig in the biggest 100 m< R bodies, for all 
numerically tested cases and for a standard dA^ oc R^^ ^dR size 
distribution. All fractions are normalized to M„,„,. 



Run 




Mhig/M„„„ 


Nominal case 


0.0678 


3770.8 


Low mass case 


0.0744 


1937.5 


<e>=0.01 


0.0318 


8515.7 


<e)=0.03 


0.0278 


3889.5 


{e)=0.2 


0.1136 


4476.2 


Hard material case 


0.0399 


2137.7 


Weak material case 


0.0840 


4112.3 




0.0639 


3655.3 


oc a"' 


0.0652 


3634.2 


oc a^^ 


0.0705 


2958.1 


dN oc R^' ^dR distribution 


0.0675 


3510.8 



inner regions than in the outer ones (see Fig.[TJ. For the inner- 
most annulus, this significant mass erosion is even observed for 
the biggest planetesimals at the upper end of the size distribu- 
tions (which get depleted by a factor ~ 2 in lO^yrs). It should 
be noted that the erosion of the ~ 0.05 mm to 1 mm grains is 
sensitive to the collisional prescription: neglecting for instance 
cratering impacts leads to a much slower evolution of this pop- 
ulation and thus a much slower flattening of the profile. 

6.2. Link between dust and planetesimals 

As described at length in the introduction, an important issue is 
the link between the observed dust population and the unseen 
bigger parent bodies. We report in Table|2]the respective masses 
of 3 representative populations: 

- the smallest R < 20yum grains, i.e., the population contain- 
ing most of the optical depth 

- all grains in the 0.1 mm to 1 cm range, i.e., where most of 
the observable "dust" mass is 

- the biggest objects in the R > 100 m range 

Surprisingly enough, the respective masses between these 3 
populations never drastically differ from their values in a stan- 
dard Dohnanyi distribution. Both M^,„/M,„„, and M^/g/M,,,,,, 
stay within a factor~3 above or below the reference values de- 
rived by integrating a dA^ oc R^^ ^dR power law. As a conse- 
quence, despite the strong wavy features of the size distribu- 
tions, the link between the total amount of observed dust and 
unseen bigger bodies can be, as a first approximation, derived 
using a simple Dohnanyi power law. 

7. Impact on the observations 

7.1. Scattered light surface brightness profiles 
7.1 .1 . Nominal case 

We consider here the two limiting cases of edge-on and 
head-on viewed systems. For sake of simphcity, we have as- 



sumed gray scattering and we display results only for the pure 
isotropic scattering case. However, other scattering phase func- 
tions have been explored, and we verify that the results pre- 
sented hereafter, in particular regarding the departure from the 
initial profiles, still hold for all explored cases. We furthermore 
assume the disc vertical scale height H varies linearly with the 
distance to the star. 

For the edge-on viewing case. Fig. [13] (left panel) shows 
that, surprisingly, the final mid-plane surface brightness (here- 
after SB) profiles only weakly vary with the parameters ex- 
plored in the different runs. For all 9 cases considered, the scat- 
tered light radial profiles approximately follow a power law in 
5Bedge(fl) °^ fl'' with -2.4 < b < -2.1. This is very differ- 
ent from what is obtained for a theoretical system where all 
bodies follow the initial S(fl) oc a" radial distribution and a 
Dohnanyi-like distribution holds over the whole size range, for 
which we get bo ^ -3. 4 (for a = -1 .5), close to the theoret- 
ical value of -3.5 (e.g. Nakano 1990f) . We shall from now on 
refer to this theoretical disc, which in fact corresponds to the 
situation at f = in our simulations, as the "static" case, with 
5Bedge(fl) 0^ fl*" and bo - a-2 (again assuming H oc a). A simi- 
lar result holds for the head-on case, for which average SB pro- 
files also strongly depart from the MMSN case (Fig. [13] right 
panel). In other words, SB profiles cannot be simply derived 
by assuming the simplest hypothesis that dust grains follow the 
same spatial distribution as larger parent bodies (for which the 
initial S(fl) oc a" profile still holds). 

Interestingly, neither can these SB profiles be derived by as- 
suming the seemingly more advanced hypothesis that all small 
(i.e. radiation pressure affected) particles have eccentric orbits 
with their periastron coinciding with the big particles distri- 
bution and their number density being derived by the classi- 
cal dA^ oc R-^^AR size distributio n. This possibili t y has been 
checked following the method of lAugereau et alJ (1200 ih and 
Thebault & Augereau I ([2005): we run a simple deterministic 
orbital integration where Q < /3 < 0.5 grains are randomly pro- 
duced from an initial parent body population (following here a 
surface density profile in a '^). The distributions of all grains 
of a given /3 are then obtained by phase mixing of their orbits 
and the total resulting surface density by weighting each con- 
tribution according to a Dohnanyi size distribution. The result- 
ing mid-plane SB profile is shown in Fig. [13] (triple dot-dashed 
line). Although it is a slight improvement over the pure "static" 
case, it is still far from all synthetic profiles obtained with our 
collisional evolution code. This means that the profile flatten- 
ing is not simply due to the geometrical spread of high-jS grains 
on eccentric orbits. It is the consequence of the more complex 
effects these movements of radiation-pressure affected grains 
have on the collision production and destruction rates of dust 
grains in the different regions of the disc. 



7.1 .2. Comparison to observations 

Different SB profiles are obtained when starting from dif- 
ferent initial radial distributions (cases explored in Sec. l5.3l l. 
However, we see that while 5'Bedge(a) profiles do vary with 
index a, the differences between the initial and final profiles 
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Fig. 13. Synthetic luminosity profiles in scattered light (isotropic scattering), for the final states of all cases explored in section 
12] (dashed lines); with the exception of those with alternative initial density profiles. All profiles have been renormalized to 1 
at lOAU. The full line represents the theoretical brightness profile should a MMSN power law in a hold, for all grain sizes, 
starting at the innermost annulus. The dashed-dotted line represents the theoretical profile obtained when assuming a MMSN 
surface density in a for all bigger grains on Keplerian orbits, and assuming that all smaller radiation-pressure affected grains, 
up to jS < 0.5, are produced by the bigger grains following a size distribution in R^^ ^ (see text for details). Left panel: disc seen 
edgeon, i.e., radial mid-plane profiles Right panel: disc seen head-on, i.e., average surface brightness. 
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Fig. 14. Normalized mid-plane surface brightness profiles in 
scattered light (isotropic scattering), for the 4 different initial 
density profiles explored in section |5] For each case, the full 
line represents the initial surface brightness profile, i.e. the pro- 
file obtained should the initial "L cc a" distribution hold for all 
sizes, and the dashed line shows the final profile at 10^ years. 



are remarkably similar, regardless of the initial S(fl) distribu- 
tion. Fig. [14] shows indeed that, for all 4 explored initial E(fl) 
distributions, the initial SB radial profiles always significantly 
flatten. The final profiles follow an approximate power law in 
5Bedge(a) ^ fl*, whose index departs from the f = case 
hy Ab = b - bo, with Ab comprised between -1 (for the 



2(fl) oc fl""-^ distribution) and -1.5 (2(fl) oc a^^ casejl. An even 
more interesting result is that, for a given system, these final 
surface brightness profiles 5Bedge ^ a* can be directly derived 
from the mass surface density distributions 2(a) oc a" through 
the following relatively simple approximate empirical law: 



5'Bedge(fl) oc fl* <-> 2(fl) oc a", with a - 2b + 3 . 



(1) 



This relation, valid for isotropic scattering, slightly de- 
pends on the anisotropic sca ttering parameter g. Assuming a 
Henyey & Greensteinl l 1941 ) phase function, we find: 



a = 2.4^7-^4.5 for | g |= 0.5 
a ^ 3.3b + 8.1 for I ? 1= 0.8 



(2) 
(3) 



A useful consequence of these relations is that they provide us 
with a tool to trace back the distribution of large parent bodies 
from the observed SB profile. It is important to point out that 
the distribution of the small grains, those dominating the opti- 
cal depth, can still be derived the "usual" way from the bright- 
ness profiles (using for example the b - a-2 relation relation 
for constant opening discs and grey scattering). The important 
result is here that recontructing the optical depth distribution is 
not equivalent to reconstructing the mass reservoir distribution. 

These results can usefully be compared to the radial lu- 
minosity profiles derived from observations. Although de- 
bris discs come in all sorts and shapes, the general ten- 
dency is that most of them have brightness profiles with a 



Let us recall that the 5Bcdgc(a) °^ a''" profile can be interpreted as 
the one the system would have in the "static" assumption (as defined in 
section [VTt . i.e., if an "equilibrium" Dohnanyi-like size distribution 
was to hold and if all particles were to follow the same spatial distri- 
bution as the largest parent-body objects (whose spatial distribution 
never significantly departs from the initial E oc a" one). 
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Fig. 15. Mean absorption cross section per unit of mass of material, k^, att - 10^ yr, for the nominal and low-mass cases (dashed 
and dotted line, resp.)- The results are compared to the situation at f = yr (R^^-^ size distribution, solid green line), while the red 
triple-dot dashed line shows a power law fit to the long-wavelength part of k^. The mean opacity obtained assuming the fit to the 
final size distribution (empirical formula given by Eq.|6]l is overplotted (blue long-dashed line). 



rather steep radial dependence in SB oc a", with typically 
-5 < b < -3.5 for edge-on discs or -4 < b < -3 for 
head-on ones (e.g lArdila e t al. 20041 Golimowski et al. 2006t 
Schneider etal] l2006: Kalas.eLalJ|2006l |2Q07|Jj. These slopes 
are significantly steeper than the typical SB^dge ^ a^'^ ob- 
tained for our nominal case with large parent bodies follow- 
ing the MMSN radial distribution in a"'^. From our parameter 
exploration, only rather extreme cases would lead to edge-on 
brightness profiles in ~ a"^^. It would require either a very 
steep S(fl) oc a^'^ surface density profile (for the unseen par- 
ent bodies) or a very high, and probably unrealistic anisotropic 
scattering parameters (| § |> 0.9). This apparent paradox be- 
tween our simulation results, which we believe are rather robust 
with respect to the flattening of the optical depth and brightness 
profile^ and observations might be understood when recalling 



that our SB^dge ^ a profile is obtained within the regions 
where a complete collisional cascade is assumed to exist, from 
micron-sized grains all the way up to big planetesimals. There 
is no obvious reason why the full radial extents of observed 
debris discs should correspond to such coUisionnaly active re- 
gions. 



we leave out of this list systems of debris "rings" with razor 
sharp outer edges, probably sculpted by gravitational perturbers, like 
Fomalhaut, HR4796 or HD 139664 

' and are moreover confirmed by preliminary simulations from 
other teams (Krivov, private communication). 



As a matter of fact, a large fraction of the luminosity ra- 
dial profiles of spatially resolved discs could correspond to re- 
gions outside the "parent body" regions of collisional activ- 
ity. For these regions outside the parent body area, preliminary 
analytical and numerical results seem indeed to show that a 
oT^ '^ slope could be a typical signature of th e presence of high 



gra i ns escaping from t heir birth region dStrubbe & Chiang 
I2OO6 : iKrivov et alj l l2006iF I. This possibility is strengthen by 
the fact that for most debris discs, the steep slopes are de- 
rived in the outer regions located at relativ ely large distances 



from the star: beyond 120 AU f or fi Pictoris ( Golimowski et al 



2004) 



lOOAU 



.2006) . 130AU for HD 15115 jArdila etal 
for HD 181 327 jSchneider et alj|2006h . ~ 40 AU for AUMic 
jKristetal.l |2005l), 55 AU for HD53143 dKalas et alJ l2006h. 
~ lOOAU for HD 32297 dSchneider et alj l2005h . Within the 
frame of the "standard" planet formation scenario, it is likely 
that these regions are beyond the limit where accretion of larg e 
planetesimals/embryos is possible (e.g iThommes et al.ll2003h . 
so that the presence of collisional cascades starting from large 
parent bodies is questionable. As a consequence, our results 
imply (within the limitations to our approach outlined in sec- 
tion [83]l that an observed SB^dge ^ oT^'^ luminosity profile is 
the signature of either: 1) an extended parent body disc with 
a sharp E oc ot'^ density decrease (Eq.[T]l or, more likely, of 



^ This outer-edge issue will be addressed in a forthcoming paper 
(Thebault & Wu, in preparation) 
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2) a region devoid of large particles beyond the main disc. 
One robust result is in any case that regions with steady col- 
lisional cascades from large parent bodies, probably cannot re- 
sult in brightness profile signatures as steep as SB edge a^^-^. 
Interestingly, for some systems where brightness profiles could 
be observationally derived in regions closer to the star, slopes 
closer to our nominal b ~ -2.2 value have been obtained. 
This is in particular true for ySPictoris where in the ~70- 
lOOAU region where most of the dust mass is believed to re- 
side, the brightness profile follows approximately SB edge °^ 
jGolimowski et al.ll2006i) . 



7.2. Thermal emission 
7.2.1 . Dust opacity 

The waviness of the size distribution is well marked for grains 
smaller than a few centimetres radius, and should have an ob- 
servational signature at far-IR, sub-mm and millimeter wave- 
lengths. The four panels of Figure[T5] show k^, the absorption 
cross section per unit mass of solid material, averaged over 
the size distribution, at four different locations in the disc. The 
curves have been obtained assuming spherical grains made of 
a silicate core and coated by water ice beyond 20 AU (see 
Sec. lB.ll for more details about the dust properties). 

In the dA^ oc R^^-^dR case, the mean opacity /c*] can be 
approximated by a power law a:" oc A^'' beyond A ~ 70- 
100 /im, with q ^ 1 for non-icy grains (a ^ asub = 20 AU), 
and q ^ 0.8 beyond cisub- These q values compare well with 
the theoretical estimates bv iDraind (I2OO6,). or the best fit values 
obtain ed for debris discs (e.g. iDent et al. I I2OOOI: iGreaves et al.l 
l2004 . Nevertheless, Fig.[l5] shows that realistic collisional 



1.0000 r 



systems do have mean opacities that strongly depart from a 
simple power law profile at long wavelengths. At represen- 
tative distances from the star (25 AU and 55 AU), the mean 
opacity ki shows a characteristic dip at A ~ 150-200 //m, 
and a bump at millimetre wavelengths for both the nominal 
and low-mass cases. At 55 AU for example, the mean opac- 
ity ratio in the SpitzerfMlPS2 and M1PS3 bands, /f70/jm/'<'i60/jm, 
amounts to 1 .7-2 times the mean opacity ratio should aR^^-^dR 
size distribution hold. Similarly, Ksio^m/Km/jm, 'fsso/jm/'cieoA/m, 
fi300/jm/*'i60/jm, are 1.7, 2.1 and 2.4, respectively, larger than 
those found for a Dohnanyi size distribution. 

In Sec. 18. II we provide an anlytical fit to the final size dis- 
tribution responsible for the waviness of the mean opacity. The 
mean opacity obtained assuming the empirical diff'erential size 
distribution given by Eq.|6] compares well to that calculated at 
our representative distance of 55 AU. 

7.2.2. Disc SED and images 

The actual impact on the disc spectral energy distribution 
(SED) is displayed in the top panel of Fig. [16] where the 
synthetic SED s have been calculated using the model of 
Augereau et al. (Il999h . The solid line on the figure represents 
the disc SED, normalized to 1 at its maximum, at f = 
(dN oc R-^ ^dR size distribution), and the bottom panel shows 
the flux ratio after lOMyr of evolution of the system. As antic- 
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Fig. 16. Disc SEDs, including both thermal emission and scat- 
tered light (which dominated over thermal emission at /I ^ 
lOjum). 



ipated, the wavy structure of the size distribution has an obser- 
vational counterpart at far-lR to millimetre wavelengths, and 
in particular a lack of emission in the 1 50-200 jum spectral 
range compared to the size distribution. The predicted 
disc colors depart from the Dohanyi case by factors that com- 
pare to the mean opacity ratios calculated above. More pre- 
cisely, the 70jum to 160 /im, 520 /im to 160 jum, 850yum to 
160 //m, and 1300 yum to 160 //m flux ratios, are 1.5-2.1, 1.3- 
1.8, 1.6-2.2 and 1.9-2.4, respectively, larger than those found 
for a Dohnanyi size distribution. 

The 1 50-200 jum spectral range clearly appears as a criti- 
cal spectral range to test the model developped in this paper. It 
requires a good sampling of the SED at long wavelengths, and 
a sufficiently precise relative photometric calibration. Several 
observational facilities working at far-IR to millimeter wave- 
lengths, will start operation in a very near future. Some of them 
are indicated in the bottom panel of Fig. [T6l to which shoul d 
be added the SCUBA-2 came ra at JCMT (.HoUand et al.ll2006l) . 
and the SOFIA observatory dBeckUnlEoOet ICaseyll2006l) . The 
PACS and SPIRE instruments onboard the Herschel space ob- 
servatory are particularly well suited to identify the dip at 
around 150-200 //m by measuri ng the exact shape of debris 
discs SEDs beyond /} ~ 70 yum (|Pilbrattll2005l: IPogUtsch etaT 



20061) . This would allow to find a direct observational signature 



of an ongoing collisional cascade in a debris disc. 

The dependence of the size distribution on the distance to 
the star, evidenced in Figs.[TO]and[TT] has direct consequences 
on the appearance of the disc, as illustrated in Fig.[T7] In the 
near-infrared, and at shorter wavelengths, light scattering by 
small (high-j6) particles dominates the disc image. The disc 
therefore shows a decreasing brightness profile with increasing 
a as discussed in Sec. 17. II In the thermal emission-dominated 
regim (mid-infrared and beyond), the disc morphology totally 
depends on the observing wavelength. At /I = 24 jum for in- 
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Fig. 17. Appearance of the high-mass disc, assumed face-on, at t 
logarithmic scale. Bottom panel: linear scale. 



10 Myr as a function of the observing wavelength. Top panel: 



stance, the disc surface brightness smoothly decreases with the 
distance from the star, while at (sub-)mm wavelengths, the disc 
shapes a ring peaked close to the outer edge of the parent-body 
disc (~ 100 AU) , a situation tha t interestingly recalls the case 
of the Vega disc jSu et al.ll2005l) . 



8. Empirical formulae for debris disc modeling 

The purpose of the present work is to numerically explore 
the collisional evolution of an extended debris disc, when tak- 
ing into account the crucial effects of impacts induced by the 
radiation-pressure affected small grains. The different results 
displayed in Secs.|5]and|2]show that, although noticeable differ- 
ences might be observed for different setups, important generic 
trends can be derived. We propose, in the following, empirical 
laws for the size distribution and collision timescales, that can 
be used for debris disc modeling as alternatives to the classical 
R^^ ^ size distribution and to the f" ,, - (rQ)"' law. 



8.1. Fit to the size distribution 

One crucial result concerns the final size distributions. For 
almost all runs, the system always quickly reaches a quasi 
steady-state, with a pronounced wavy distribution which 
strongly departs from a standard "equilibrium" distribution in 
dA^ oc R^^ ^dR, or any simple power law in R''dR for that mat- 
ter. As clearly appears in Figs.[T]&[3] the waviness varies with 
location in the system, it is less pronounced close to the inner 
edge, since it is mostly due to collisions due to high velocity 
outward moving small grains. However, if ones considers the 
average distribution integrated over the whole disc, we have 
seen that its profile only weakly depends on parameters such as 
the system's total mass, it's dynamical excitation or the value 
of the radiation pressure cut-off size 7?pr. For the latter case, 
what is observed is mostly an offset of the wavy-distribution, 
which retains its global shape and main characteristics. As for 
the total initial mass, it does not crucially affects the final shape 
of the size distribution as long as collision lifetimes of dust 
grains are shorter than the system's age (see Sec. l4.2l i. The fi- 
nal size distribution is even relatively unaffected by the profile 
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Fig. 18. Dashed lines: Final size distribution profiles, averaged 
over the whole system, for all numerically tested cases, except 
the 2 "weak" and "hard" material cases with different colli- 
sional prescriptions. All profiles' jc-axis have been renormal- 
ized to units of /?pr, and all >'-axis to the value of the first wavy 
"peak" at /? ^ 1 .5Rpr. Solid line: reference profile derived from 
our empirical fit given by Equ.|4]and integrated over one loga- 
rithmic size interval: AM(R) = G(/?) (/^/I.S/^pr)""*' AM(i.gRp^). 



of the initial mass distribution (exponent of the S(fl) profile). 
The only cases for which a major modification of the size dis- 
tribution is observed are the "very weak" and the "very hard" 
material cases. Apart from these 2 exceptions, for all other 10 
tested setups we obtain very similar features: a strong depletion 
of R ^ /?PR grains, a peak for R ^ I .5/?pr followed by a deep 
depletion of objects in the 107?pr < R ^ 507? pr range. The sim- 
ilarities between all profiles are even more striking when they 
are renormaUzed by their value atR = /?pr (Fig. [18). As can be 
clearly seen, variations are very limited for R $ 100/?pr. For 
this size range it seems thus reasonable to consider that, as a 
first approximation, the size distribution obtained in our nom- 
inal case is a relatively good standard for spatially extended 
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Fig. 19. Dashed lines: Collisional lifetimes at aq - 55AU, 
for all tested cases except the weak and hard material runs, 
when renormalized by the reference timescale fcoii-ief = 
?"^jj(ao) (Q(ao)'!"(flo)) '■ Solid line: profile derived from our em- 
pirical fit (Equs.|7]and[8]i 

systems. We were able to derive an empirical fit for this revised 
size distribution, valid in the R ^ 100/?pr range. When written 
in terms of the differential mass distribution, it reads: 



dM oc G{R)R-^ ^'^dR 
with 

logio(GW)= ^ 
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(5) 



This new relation proves to be a reasonably good fit to almost 
all profiles in the R ^ 100/?pr range (Fig.fTSll. In terms of the 
differential size distribution dN(R), this translates into 



dA^ oc GiR)R'^-^'^dR , 



for -RpR <R<,lOQRpR. 



(6) 



Beyond 100/?pr, stronger divergences between different runs 
are observed. However, as a rough first order approximation, 
the differential size distribution can approximately be extrap- 
olated by a R^^ ^ power law. The R^^-^ extrapolation has been 
used to calculate the mean opacity represented by a blue long- 
dashed line in Fig.fTSl 

8.2. Fit to the collisional particle lifetime 



As shown in Sec. 14. 31 colUsional lifetimes strongly vary with 
particle sizes: they increase very rapidly when R gets close 
to the blow-out limit Rpr, reach a sharp minimum around 
R ^ lOiJpR, increase sharply again between 107?pr and about 
100/? PR and then continue to increase much more slowly with 
increasing sizes (see FigHJi. We have also shown that a direct 
consequence of this result is that collisional lifetimes cannot be 
directly derived from the optical depth through the simplified 
formula fi^nia) = (tQ)"'. There are several reasons why this 
formula cannot hold here: 



- the f[?(,,i(fl) = (TO.y^ formula implicitly considers impacts 
between objects of equal sizes, thus neglecting the broad 
size spectra of all possible impactors on a given target, 

- it also implicitly assumes that all impacts are fully destruc- 
tive, i.e., that the collision timescale is equal to the colli- 
sional lifetime. This neglects all cratering impacts, whose 
role is crucial for the considered problem (see FigEl), 

- even more important: this formula neglects all effects due 
to the specific dynamics of the smallest grains affected by 
radiation pressure, 

- last but not least: at any given distance ao from the star, 
it neglects all collisions due to objects coming the inner 
a < ao regions, and it has been shown (FiglJI that these 
collisions are crucial for the evolution of dust grains. 

As can be seen for example in Figs.|4^,b and |6^,b, colli- 
sional timescales significantly vary for different initial condi- 
tions, in particular total initial mass and dynamical excitation 
of the system. However, the profiles of the tcoiiia,R) curves 
are relatively similar In order to visualize these similarities 
more clearly, all tcoiiia,R) curves have been renormalized by 
the reference timescale t^^i^ia) - (tQ)"' (FigfT9]l. In a similar 
fashion as for the size distributions, we see that all normal- 
ized tco\i{a,R) profiles remain relatively close to the nominal 
cas^. This opens the possibility for deriving an empirical fit to 
fcoii(a, R) as a function of a and t: 



(4) fcoll(fl,-R) = feoll(«) 



R_y^ iRf 



for R<R2 



with/?! = 1.2/;pR andRo = IOO/^pr, and 



R\' 



tcou(a,R)^t%{a)\j^\ forR>R2 



(7) 



(8) 



8.3. Approximations and limitations 

Let us state again that these relations should be taken with care. 
An important general remark is to again stress that they have 
been derived for extended coUisionally active regions, i.e., re- 
gions with steady collisional cascades starting from large reser- 
voirs of big unseen parent bodies. These regions might not ac- 
count for all the observed radial extent of debris discs: some ob- 
served regions are probably coUisionally inactive areas where 
only small high-yS grains, produced in parent body regions fur- 
ther inside, are present (see discussion in Sec. 17.1.21 ). 

Moreover, within the frame of our numerical approach it 
is important to point that these fits are valid for our nomi- 
nal collision outcome prescription, and that significant vari- 
ations should be expected for harder or weaker material pre- 
scriptions (Fig.|9]l. It should also be noted that in a "real" disc, 
all individual particles are not completely identical: they would 
have slightly different material compositions, porosities, differ 
in presence or absence of microcracks, etc... This might alter 



^ The significant difi"erences observed between the dynamically ex- 
cited and dynamically cold cases (see Fig[6l( are partially erased after 
renormalization by fj^^j/a)- Indeed, as seen in Fig[5] systems with low 
(e) are globally depleted in R ^ 0.1mm grains and have thus lower 
optical depth (since t is mostly contained in the smallest particles) 
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the size distribution profile, probably damping the waviness 
described in Eq.|6]to some extent, but such sophisticated ef- 
fects are difficult to take into account with a particle-in-a-box 
code. Another important point is the fact that the smallest par- 
ticles considered here are just below /?pr, so that only 2 size 
"bins" correspond to unbound so-called "y6-meteoroids". We 
nevertheless performed a few test runs with additional small- 
size bins, and observed no drastic change in the final pro- 
files. However, for more massive discs, taking into account the 
role of j g-me teoroids, as was done in the pioneering work of 
Krivov et al. might be crucial. For such high-mass sys- 

tems, extremely efficient collisional "avalanches" chain reac- 
tions trig gered by /6-meteoroids could possibly play a signifi- 
cant role (IGrigorieva et al ]|2007h . The contribution of unbound 
grains could also be important for interpretation of observa- 
tions particular ly sensitive to smaller particles, e.g. polarimetry 
jKrivova et aLirioOOl) . 

We do however believe that, regardless of their exact level 
of accuracy, the present empirical fits are in any case a more 
reliable fit to "real" size distributions than any simple dA^ oc 
R'^dR power law (be itq- -3.5 or not) extrapolation. 

9. Summary and conclusions 

We elaborate in this paper a model able to follow the collisional 
evolution of extended debris discs ov er a 10 Myr span. We c on- 
firm the previous results obtained bv iThebault et al.l (l2003b for 
a narrow, isolated annulus, that the classical dA^ oc R^^-^dR 
Dohnanyi size distribution cannot hold in realistic collisional 
discs. Rather, a wavy size distribution develops in the whole 
system, amplified by the particular dynamics of the radiation 
pressure affected grains (high-yS particles). 

The model builds on the classical particle-in-a-box tech- 
nique, and allows a detailed exploration of the various parame- 
ters that impact the disc evolution. Such a quantitative numeri- 
cal exploration had not been undertaken so far, at least not when 
following the size distribution evolution over a range encom- 
passing all objects from the fim to the biggest parent bodies in 
the 50 km rang^ We chose therefore not to focus on one given 
observed debris disc but to consider a fiducial nominal system, 
making the most reasonable (or maybe least unreasonable) as- 
sumptions, in order to clearly identify and quantify the complex 
mechanisms at play, and derive general behaviours without bi- 
ases by non-generic artifacts. However, in order to check the 
robustness of our results, several key free parameters have been 
explored. Our main results can be summarized as follows: 

1 . A wavy size distribution, strongly departing from a dN oc 
R^^-^dR power law, is a common feature of collisional de- 
bris discs. 

2. The wavy pattern includes an overdensity of grains with 
radius about twice the blow-out grain size RpR, and a strong 
depletion of the 10 - 507?pr particles. 

3. The waviness weakly depends on the disc mass, initial sur- 
face density profile, mean disc dynamical excitation, stellar 



with the notabl e exception of the v ery innovative and promising 
kinetic approach of lKrivov et al.l ( l2006h . but so far considering a very 
simplified model of collision outcomes 



properties, but is affected by the collision outcome prescrip- 
tion, especially the resistance of objects to collisions. 

4. In extended discs the evolutions of different regions of the 
systems are strongly interconnected: the waviness is ampli- 
fied by high-y6 bound particles (grains strongly affected by 
pressure forces), which have large radial excursions within 
the system and can impact, at very high velocities, larger 
objects far outside the region where they were initially pro- 
duced. 

5. Surprisingly, the global dust to planetesimal mass ratio is, 
to a first order, not strongly affected by the size distribution 
waviness. 

6. Collisional lifetimes strongly differ from the usual (tQ) ' 
approximation in realistic collisional systems. 

7. The optical depth and the scattered light flux are domi- 
nated by a very narrow range of so-called a-meteoroids, 
i.e., bound objects just above the blow-out cutoff size. 

8. Spatial distributions are also affected. The radial distribu- 
tions of grains of different sizes might significantly diverge 
from one another. More generally, there is a major dis- 
crepancy between the radial distribution of particles in the 
dust-size range, i.e. those accessible to observations, and 
the largest undetectable objects that make up most of the 
system's mass. The distribution of small grains, and thus of 
the disc's optical depth, is significantly flatter than that of 
the big parent bodies. 

9. This flattening of the small grains radial distribution trans- 
lates into a flattening of surface brightness profiles in scat- 
tered light in the regions where the big parent bodies reside. 
For a disc having an initial MMSN surface density profile 
the equilibrium scattered light surface brightness profile is 
roughly in SBgdge °^ cfi, with -2.3 < b < -2 instead of the 
standard b ^ -3.5 value. 

10. These radial slopes are less steep than those observed for 
the vast majority of debris discs. This apparent paradox 
could be explained by the fact that for most systems, ra- 
dial brightness profiles are observed in regions beyond the 
outer edge of the main "parent body" disc. In these regions, 
no collisional cascades take place and only small high-/? 
grains, produced further inside and pushed on eccentric or- 
bits by pressure forces, are observed. 

1 1 . The waviness of the size distribution translates into wavy 
dust opacities and SEDs at far-IR and (sub-)millimeter 
wavelengths, which could be observable signatures of the 
collisional activity in debris discs. 

12. We derive an empirical formula for the differential size dis- 
tribution (Eq.|6]l which fits reasonably well the numerically 
obtained results. Although this approximate fit should be 
taken with care because of the unavoidable limitations of 
our numerical code, future models aiming at reproducing 
multi-wavelength observations might use this formula as 
an alternative to simplified dN oc R^dR power laws. 

13. Similarly, we propose an empirical formula for the colli- 
sional lifetime of the particles (Eq.|2] & [8]l that might be 
used to interpret data. 

This paper provides the basis for future debris discs mod- 
eling of individual cases such as Vega, for which both resolved 
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data and numerous photometric measurements are available. 
But overall, the waviness of the size distribution is becoming a 
well established feature that cannot be ignored in future SED 
analysis, and the empirical size distribution given by Eq.|6]is 
provided for this purpose. We in addition stress that a wealth of 
future facilities working at far-lR and (sub-)millimeter wave- 
lengths (Herschel, SOFIA, SCUBA-2, ALMA) will soon of- 
fer the opportunity to test the model developed in this paper, 
providing a direct observational hint for an ongoing collisional 
cascade in a debris disc. 
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Appendix A: Evolution equation 

The present muhi-annulus code is based on the single-annulus 
algorithm developed in TAB03 and described at length in this 
paper We shall thus only recall here its main characteristics 
before describing in more details the enhancement performed 
for the present version. 

We first spatially divide the system into concentric an- 
nuli. Within each annulus, we follow the classical particle-in- 
a-box approach in which the particle population is divided into 
n boxes each standing for a given particle size R^. In each an- 
nulus ia > 2 (all except the innermost one), additional bins are 
included which account for the small grains originating from 
ia' < ia annuli and placed by radiation pressure on highly ec- 
centric or unbound orbits crossing the ia annulus. We arbitrar- 
ily set the limit for grains sizes for which additional bins are 
considered by the criteria > /3\im - 0.05. For one given parti- 
cle size Ri in the ia annulus, there are thus 1 + «;,(,) correspond- 
ing bins, where < nb{i) < ia is the number of possible source 
annuli ia' < ia for all "foreign born" Ri grains. To describe 
the number of particles of one given grain population within 
one given annulus ia, we use the terminology Niaj^ia', where 
1 < ! < n is the size bin index and ia' < ia the source annulus 
where the grain population has been produced {ia' = ia for all 
particles with /3i < jSum)- 

At each time step, the change in the number Niajja' is given 
by the collisional evolution equation displayed in Equ.l of 
TAB03. For the locally produced small grains affected by ra- 
diation pressure > /Sum), an additional term is introduced 
which reads 



foi 



out(iaJ,ia) ^iaJJa Jin{iaJ,ia) ^"^ia+ljja 



+ fin 



(A.l) 



where fout(ia,i,ia) is the fraction of N/ajja particles leaving the ia 
for the ia + 1 annulus during dt and fin(ia,i,ia) the fraction of par- 
ticles re-entering the ia annulus after completing one full orbit 
{fin{ia,k,dr) - for grains on unbound orbits). Note that these 
re-entering particles necessarily come from the neighbouring 
ia+l annulus, where they were members of the Nia+ijja bin. All 
fout{ia,i,ia) and fi„(ia,ija) ratcs are derived from separate determin- 
istic numerical simulations following the dynamical behaviour 
of 10000 test particles with = y6, released on randomly dis- 
tributed orbits with e = 2/ = (e)o and a^in = IQAU < a < 
amax = 120A;7. 

For the y6, > 0.05 grains which have been originally pro- 
duced in an inner ia' < ia annuli, the additional evolution term 
due to inter-annuli exchanges reads 



ia.i,ia' 



(A.2) 



where gout+mM') and gou,-(ia,i,ia') are the fraction of outgoing 
(to the ia + 1 and ia - 1 annuli respectively) particles and 
ginHiajM') and gin-{ia,iM) the fraction of incoming (from the 
ia + 1 and ia - 1 annuli) particles. The terms gin-(ia,ija') and 
gout+{ia,ua') Correspond to particles produced in the ia' annu- 
lus on their way out towards their apoastron (or infinity for 



unbound orbits) and the terms gi„+(iaj,ia') and goui-(ia,Ua') cor- 
respond to particles having already reached their apoastron and 
on their way back to the ia' annulus (these terms are equal to 
zero for unbound orbits). All 4 parameters are estimated with 
the same type of numerical simulations as those used for deriv- 



ing/. 



ouHiaJJa) 



and fi. 



m{ia,i.ia) ■ 



As already mentioned, the dynamical state of the system is 
fixed and does not evolve with time. To estimate the average 
encounter velocities, we divide all possible target-impactor en- 
counters into two types: 1) those involving two < 0.05 (i.e. 
not significantly affected by radiation pressure) particles, and 
2) those where at least one of the involved bodies is on a radi- 
ation pressure modified orbit (y6, > 0.05) For type 1) impacts 
within one given annulus ia at a distance r,v, from the star, the 
enco unter velocity is sim ply given by the classical expression 
(e.g. lLissauer & Stewarl l993; Thebault et alj|2003b : 



ia;jja 



1/2 



(A3) 



where {vkep(ia)) is the average Keplerian velocity at distance r,fl, 
and (e) and (/) are the average orbital parameters imposed as 
initial conditions. As described in section 2, we take here (e)o - 
0.1 = 2 (/)o. For type 2) impacts the average impacting speed 
are numerically estimated in specific determinisitic numerical 



Appendix B: Collision outcomes 

We follow here the classical approach where collision out- 
comes are divided into 2 regimes, depending on the ratio be- 
tween the specific impact energy per target mass unit Qim^ - 
Eco\/Mt, where E^oi - MpMtAv'^/2(Mp + Mt), and the critical 
specific energy Q*: catastrophic fragmentation if Qimp > Q* 
and crate ring if Qimp < Q*- 



B.1. Critical specific energy 

Q* is a function of the target's radius R,, this dependence being 
usually expressed as th e combination of two power laws (e.g. 
Benz&Asphau3ll999l) : 



^ " \ 1cm/ \ 1cm/ 



(B.l) 



The first term on the right hand side corresponds to the strength 
regime, valid for small sizes, where Q* slowly decreases with 
size, while the second term has a positive index correspond- 
ing to the gravitational binding regime. Values for Qqs, B, 
bs and bg depend on the physical composition of the ob- 
jects and are derived from laboratory experiments or numerical 
models (e.g. iHousen & H ol saDplel[l990l: lOavis & Rvanl | l99C : 
Holsappldl 19941; IPaolicchi et al.l ll996l:lBenz & Asphaug|ll999 : 



Arakawall 19991) 



A first important issue is then which chemical composi- 
tion to assume for the objects: proportion of ices and silicates, 
porosity, etc... For our nominal case of a A5V star, we have as- 
sumed simple mixtures of silicates, water ice, and vacuum (to 
mimic porosity). The ice subHmation distance is a function of 
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Fig.B.l. Position of the snow line (T =^ 120 K) for 
grains with silicate-rich cores around a jS Pictoris-like star. A 
NextGen synthetic stell ar atmosphere spectrum for a A5V star 
jHauschildt et al. I I1999I) las been used to compute the equi- 
librium temperature of the grains. The silicates and water ice 
optica l constants are from IPraind (l2003l) and iLi & GreenbergI 
(119981) . respectively. The mean grain optical index was ob- 
tained using the Maxwell-Garnett mixing rule, and the Mie 
theory used to compute the grains absorption/emission efficien- 



the grain size, as shown in Fig. lB.ll for two different vacuum 
volume fractions (0% to simulate non-porous grains, and 90% 
for highly porous grains). For grains larger than a few jim, the 
sublimation distance oscillates about 20+5 AU, with the largest 
distances reached for the smallest grain sizes considered in this 
paper Given the spatial resolution of our simulations, we adopt 
a single, average sublimation distance for water ice of ~ 20 AU. 

But the problems do not stop here since, even for similar 
materials, the different Q* prescriptions available in the liter- 
ature do often sig nificantly diverge from o ne another (see for 
example Fig. 8 in iBenz & Asphauj 1999). For icy bodies in 
particular, critical energy estimates might differ by up to two 
orders of magnitude dependi ng on the studies (see for example 
the discussion in section 4 of Burchell et al ] l2005h . It is not the 
purpose of the present work to address this very difficult issue. 

We shall here consider a "nominal" case, assuming 
that for sihcates the cr itical energy is the one derived by 
Benz & Asphau g|(ll99^ for impacts at 3km.s i.e., Qos = 



3.5 X lO^erg.g"', B = 0.3 er g.cmlg' ^, b, = -0.38 and 



bg - 1.36. For ices, we follow i Krivov et al. ( 2006h and take 
2ice{R) ~ ^/^2sii(R)' which is an approxim ate intermediate 
positio n between the hydro-code results of ,Benz & Asphaug 
1999h . who found that ice could be almost as resistant as sil- 



icates, and most experiment results in which ic es proved to be 



more than one order of magnitude weaker (see iBurchell et al 



2005L and reference therein). Because the smaller radiation- 
pressure affected grains might impact objects in regions differ- 
ent from the ones where they have been produced, we have 
also to take into account impacts where one colliding body 
is icy and the other one rocky. For these heterogeneous colli- 
sions, we assume a simple prescription with Q*^^_^^^ = 1/22*,^. 
and 2sii-ice ~ ^Q*siv roughly taking into account the fact that 
impacts by hard (resp. weak) projectiles on weak (resp. hard) 
targets are more (resp. less) destructive than impacts between 
bodies of same material. Furthermore, since impacting veloci- 



ties do significantly vary within the disc, we take into account 
the we ak Q* dependence on Av found by lHousen & Holsapple 
( 1990l) and assume 



<2(A.o = e; 



(Sim.j-') 



Av 



3km. s 



0.35 



(B.2) 



Note that our critical energy prescription g i ves Q * values 
significantly higher than those of iKrivov et al. ( 20061) . These 
authors assumed 2o.s = 3 x 10''erg.g~^ (at R=lcm) for silicates, 
which is significantly belo w most Qnt estimates a vailable in 
the literature (see F ig. 8 of |Benz & Asphaug|[l999l) . with the 
exception of that of Durda et alj ( 1998 ) obtained from fitting 
the observed size-d istribution of asteroids. The prescription of 
Krivov et al. I (l2006h will however be tested as a "weak mate- 
rial" run (see section l54b . 

B.2. catastrophic fragmentation 

If 6imp > Q* catastrophic fragmentation occurs: the target is 
shattered and produces a population of fragments where the 
biggest one has a mass M\f < 0.5 M,. The value of M\f as well 
as the size distribution of the produced fragments is computed 
following the procedure described at length in section 2.4 of 
TAB03. 

We would like to point out that our model departs from 
the often assumed simplifying assumption that fragment size 
distribution follows a power law in d N oc R^^ ^dR (e.g. 
Augereau et"al]|200ll: IkjIvov et alj|2006i) . Such a power law 
is indeed in principle the equilibrium value reached after suffi- 
cient mutual collisions, but not the one for fragments produced 
after one given impact. Furthermore, we consider here a broken 
power law, with 2 different indexes, corresponding to a change 
of slope for the size distributions of the smallest fragments, 
a feature w hich is supported by experimental and theoretical 
studies (e.g. bavis & Ryanlll990l:lTanga et alJll999l) . 



6.3. cratering 

If Gimp < Q*, the target is preserved but eroded by a mass 
Mcia. In most published collision-evolution models, Mem is 
directly proportional to E^o] through a constant coefficient 
g, oft en called excavation coefficient (e.g. Greenberg et al.l 
1 19781: LStoffleretanil975l: IPetit & Farinellalll993b or defined 
as a = i/Qc, where Qc is the "crus hing energy" (e.g. 
Wetherill & Stewart! 1 1993UKenvon & Luu| [l999). Values of a 
typically range between 10^'s^cm"^ for hard material and 
4 X 10-**- 



^ for weakly b onded sand (e.g. 'Greenberg et al.l 

19781 iDobrovolskis & Burrisi ri984: Petit & Farinella 1993). In 
TAB03 we also followed this prescription, with an intermedi- 
ate a = 10"**s^cm"^ value. However, the Mem = ff^coi relation 
is in reality a simplification of the more general 



Me, 



-'EL 



(B.3) 



dependence, with index j slig htly greater than 1 (e.g. Gaultl 



19731: iKoschnv & GriinI l2001i) . 



Mcra - aEeoi relation has been derived by iMarcus 



Historically, the simphfied 
19691) 
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who extrapolated experimental results, obtained mainly by 
Gault et alj ( Il962h on small projecti les, to the much larger sizes 



with 



considered in his study (see p. 77 of Marcusl[l969h . It has been 
later assumed bv lOreenberg et alJ ( 119781) in their milestone nu- 
merical study of planetesimal accretion and in most statisti- 
cal collisional evolution models ever since. But one should be 
aware that this relation is in principle only valid over a limited 
range of object sizes and velocities, typically lOm^ R ^ 1km 
and Ay ^ 3-5 km.s"', and that a study considering size ranges 
spanning over several orders of magnitudes should assume the 
'real" dependence in a'E^ .. One of the m ost accurate Mg 



perscr iption is probably given in Equ.7 of i Koschny & GriinI 
(1200 ih . which is an empirical fit of experimental results ob- 
tained by these authors for mixed ice/basalt bodies as well as 
by several other studies for pure silicate or pure ice objects and 
reads (with the present formalism); 



logiK) ^ logiK,,) - 



log{Q.2) - logiQ.Ql) 
where 



(log(K,,) - log(Ku)) 



Q* 



(B.8) 



(B.9) 



In a similar way to what was assumed for the Q* parameter, 

we assume that Mcra(ice) = 5 Mci-a(sil) as well as Mcra(ice-sil) - 
2Mcra(ice) and Mcra(sil-ice) - 1 /2 Mcra(sil)- 

The excavated mass Mem is then redistributed into frag- 
ments following a single-index size distribution power law (see 
section 2.5 of TAB03). 



sil 



M„a = Vice I TT^ I 2-V£ 



(B.4) 



where fsi\ is the proportion of silicates in the target, Vsn = 
lO^'^cgs, Vice = 6.69 X lO^^cgs, and y = 1.23. Note that this 
formula gives substantially lower excavated masses for small 
(< 1cm) targets than those derived with the M^a - 10"'*£'coi 
relation taken in TAB03. 

Nevertheless, this formula is only valid in the small- 
scale impact regime, where M^a << Mt, corresponding ba- 
sically to a grain-hitting-a-wall case. For larger craters, ef- 
fects of cratering in finite spheres have to be taken into account 
(IHolsappldl 1 994t) . This raises the more general issue of "con- 



necting" the cratering prescription to the fragmentation one. 
Some collision-evolution models assume an abrupt fragmenta- 
tion/cratering transition, where the maximum possible value of 
Mg-d/Mt just below the fragmentation t hresh old is ^ 0.1 (e.g. 
Petit &Farinellal 1 19931: iThebault et aP l2003h . thus implicitly 
leading to a sharp drop from (1 -Mif) - 0.5M, to M^a - 0. IM,. 
Experiments seem however to show that there is no sharp tran- 
sition around the M |f = 0.5M, value dPavis & Ryan Il990t 
Housen et al. 199 l l). so that the transition between the frag- 
mentation and cratering regimes should be more or less pro- 
gressive. Such a smooth transition is assumed for our present 
model, where we consider 3 cases. For small-scale craters, we 
take: 



M„a = a'£j„i if gimp < 0.01 e* 



(B.5) 



with a' given by Equ lB.4l For the large- scale regime jus t belo\y 
the fragmentation threshold, we follow IWvatt & Dent I (120021) 
and assume: 



M„a = 0.5M, 



imp 



Q* 



if 0.2G* < a„p < Q* 



(B.6) 



which is in agreement with th e experiment results displayed in 
Fig. 5 of iHousen et alj (1199 11) . Between these two modes, we 
assume a smooth transition given by: 



M,,,^KE,,i if 0.01 e* < anip < 0.22* 



(B.7) 
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